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ABSTRACT 

Large  matrix  storage  constitutes  a  limitation  on  the 
applicability  of  most  numerical  techniques  including  the 
Finite  Element  Method,  when  very  accurate  results  are  required. 
This  is  particularly  true  when  dealing  with  Boundary  Value 
Problems.   In  order  to  surpass  this  difficulty  a  new  method 
to  solve  these  problems  has  been  devised  which  does  not  re- 
quire matrix  storage  while  still  providing  the  possibility  of 
accuracy  improvement. 

Although  restricted  to  one-dimensional,  linear  differen- 
tial  equations  of  the  form  Y  (x)  =  f (x)  this  new  approximating 
technique  gives  acceptable  results.  The  method  will  perform 
equally  well  for  problems  with  exact  or  non-exact  integrable 
forcing  functions,  continuous  or  discontinuous,  or  functions 
existing  only  as  a  set  of  values  at  discrete  points. 
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I.   INTRODUCTION 

For  many  engineering  problems  it  is  not  always  possible 
to  find  an  exact  solution.   An  exact  solution  is  an  analyti- 
cal mathematical  expression  that  gives  the  value  of  the  unknown 
at  any  point  within  a  previously  specified  range.   For  problems 
involving  complex  material  properties  and  boundary  conditions, 
different  numerical  techniques  have  been  developed  to  approxi- 
mate the  exact  solution  to  a  more  or  less  acceptable  degree 
of  accuracy.   One  of  these  techniques  is  the  Finite  Element 
Method,  where  a  given  problem  is  to  be  discretized,  that  is  to 
say,  approximate  solutions  are  to  be  found  at  discrete  points 
in  the  body.   The  way  one  accomplishes  this,  is  by  subdividing 
the  whole  body  into  finite  elements.   The  solution  is  then 
formulated  for  each  small  unit  and  combined  to  obtain  the 
solution  of  the  whole  system.   Obviously  the  greater  the  number 
of  elements,  the  better  the  accuracy. 

However,  despite  the  fact  that  high-speed  digital  computers 
have  enabled  engineers  to  successfully  apply  these  numerical 
techniques,  we  still  face  the  problem  of  large  matrix  storage. 
Complex  problems  requiring  very  accurate  approximate  solutions 
lead  to  large  matrices  and  consequently  larger  computer  memory 
will  be  required.   This  fact  constitutes  a  limitation  on  the 
applicability  of  these  numerical  techniques. 


The  method  we  are  dealing  with  in  this  work  intends  to 
solve  this  difficulty  by  reducing  a  B.V.P.  to  an  I.V.P.   Unlike 
the  "shooting"  method,  where  basically  the  same  idea  is  used 
together  with  an  iterative  scheme  to  achieve  a  solution,  the 
method  developed  here  achieves  a  solution  without  iteration. 

Although  restricted  to  a  particular  type  of  linear  differ- 
ential equations  and  only  to  one-dimensional  problems,  as  it 
will  be  shown,  this  method  can  be  applied  to  many  problems 
where  no  exact  integrable  forcing  functions  occur,  or,  even 
more,  the  function  exists  only  as  a  set  of  values  at  specific 
points . 

Basically  what  the  method  does  is  approximate  the  exact 
solution  by  a  set  of  linear  functions,  each  of  them  applying 
in  a  very  small  interval.   However,  these  functions  are  not 
independent  from  each  other.   We  use  the  previous  one  to  find 
the  next,  until  we  eventually  reach  the  other  end  of  the  range 
in  consideration.   Depending  on  the  order  of  the  differential 
equation,  we  will  find  as  many  sets  of  linear  functions  as 
required  by  its  order,  each  set  applying  to  a  specific  inte- 
gration, and  as  before,  every  set  depends  on  the  previous  one. 
The  way  in  which  we  will  do  this  will  allow  us  to  correlate 
all  possible  boundary  conditions  in  an  explicit  manner  which 
transforms  boundary  value  problems  into  initial  value  problems 
and  by  doing  that,  get  to  the  result. 

The  method  is  indeed  an  approximation  and  as  such,  it  is 
subject  to  error.   We  can  say,  however,  that  by  decreasing  the 


step  size  we  get  better  accuracy,  provided  we  carry  enough 
significant  digits  to  take  care  of  round-off  errors.   With 
respect  to  the  later  a  good  programming  technique  is  required. 
Finally  we  shall  apply  the  method  to  several  practical  applica- 
tions, specifically  to  beam  problems  where  a  fourth  order 
differential  equation  occurs,  and  a  variety  of  boundary 
conditions  can  be  given. 


II .   BACKGROUND 

A.   GENERAL 

In  this  section  we  are  going  to  introduce  the  basic  ideas 
of  the  method  we  will  be  dealing  with.   Let 

Y'(x)   =   f(x)  (2.1) 

where  it  is  understood  that  the  function  f (x)  may  or  may  not 
be  exactly  integrable,  but  we  can  integrate  it  numerically. 
From  the  above  relationship  we  have: 

Y(x)   =  /f(x)dx  +  YQ  (2.2) 

Here  the  constant  of  integration  Y   is  assumed  to  be  zero 
for  the  moment.   In  the  most  general  case,  it  is  clear  that 
there  is  no  way  to  know  what  the  function  Y(x)  is,  since  we 
may  be  dealing  with  non-integrable  functions.   However,  we 
can  approximate  the  function  Y(x)  by  a  linear  function  y(x), 
provided  the  interval  in  which  this  approximation  applies  is 
sufficiently  small.   So  we  take 

Y(x)   =y(x)   =mx+c  (2.3) 

Now,  by  following  the  general  procedure  of  integration,  and 
from  (2.2) ,  we  have: 


b       b 
mx+c|    =    /   f ( x ) dx  (2.4) 

a     a 


where  (a,b)  denote  the  limits  of  the  interval  under  considera- 
tion.  After  simplification  it  is  found  that: 

b 
/   f(x)dx 

m   =   - — r- (2.5) 

b  -  a 


that  is  to  say,  the  slope  of  the  approximating  line  given  by 
(2.3)  can  be  explicitly  determined.   Actually,  as  we  will  see 
later,  it  can  be  said  that  the  "general"  approximate  solution 
to  (2.1)  has  been  found  within  the  interval  from  a  to  b .   In 
order  to  determine  the  "particular"  solution,  in  other  words 
to  find  the  intercept  c  of  (2.3) ,  an  auxiliary  condition  must 
be  imposed.   Let 


Y(a)   =  Y 
o 


then: 


Y(a)   =   Y    =   ma  +  c 


or 


y(x)   =   m(x-a)  +  Y0       a  <  x  _.  b  (2.6) 
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where  m  is  given  by  (2.4) .   Now  we  have  fully  determined  the 
approximate  solution  inside  the  given  interval.   Note  that 
(2.6)  will  give  results  at  points  a  and  b,  as  acourate  as  the 
numerical  integration  performed  on  (2.5). 

Fig.  2.1  shows  the  whole  process  so  far.   The  upper  curve 
represents  f(x),  which  is  given.   The  lower  curve  is  Y(x), 
the  integral  of  f(x);  however  within  the  interval  a  to  b,  Y(x) 
is  approximated  by  a  line.   Point  Y   is  the  given  auxiliary 
condition  and  point  Y,  can  be  determined  exactly  from  (2.6), 
by  letting  x  =  b. 

The  next  step  now  becomes  evident,  since  point  Y,  can  be 
determined.   We  are  now  in  the  same  position  as  before,  so 
all  we  need  to  do  is  repeat  the  process  over  again.   However, 
this  time  with  a  new  auxiliary  condition;  the  last  one  we 
have  just  found,  and  over  the  next  interval.   We  keep  going 
this  way  until  we  reach  the  other  extreme  of  the  range  where 
the  differential  equation  applies.   See  Fig.  2.2. 

In  summary,  we  can  say  that  every  step  we  take  we  are 
solving  a  "new"  differential  equation  by  approximating  lines 
which  gives  us  true  values  at  the  points  of  intersection. 
The  fact  that  we  have  a  "new"  differential  equation  at  every 
single  step  allows  us  to  deal  with  discontinuous  functions  and/ 
or  functions  existing  only  as  a  set  of  values. 

B.   THE  FIRST  INTEGRATION 

At  this  point  let  us  introduce  a  new  parameter: 

h   =   b-a  (2.7) 
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f(x) 


Yw=fcx) 


y(x)  =  f[(xU: 


b  % 


FIGURE     2.1 
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*1  *3  *n-t 


FIGURE    2.2. 
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where  h  is  assumed  to  be  small  and  constant  throughout  the 
analysis.   We  are  now  in  a  position  to  determine  a  general 
relationship  which  will  allow  us  to  find  the  value  of  Y(x)  at 
discrete  points,  namely  at  the  intersection  of  the  straight 
lines.   From  (2.6)  we  have: 


y(b)   =   m(b-a)  +  Y    =   mh  +  Y 
J  o  o 


In  general, 


where 


but : 


Y.  ,   =   m.h  +  Y.  ,  0<i<n         (2.8) 

l-l       l      l-l       —   — 


x . 

l 


/    f(x)dx 


x .  , 

m.   =   -i^± r- (2.9) 

l  h 


Y.  ,   =   m.  ,h  +  Y.  0 
l-l       l-l      i-2 


Y   „   =   m.  „h  +  Y.  -, 
i-2       1-2      1-3 


Y,   =   m,h  +  Y 
1      1     o 


Replacing  these  last  relationships  in  (2.8)  and  after  factor- 
ing h,  we  have: 
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Yn      =       (m,    +   iru    +    .  .  .    +   m    )h    +   Y 
n  12  no 


or 


n 


rn      =      h      I      m      +   Y  (2.10) 

i=l      1  ° 


Similarly, 


m 


h 
/      f(x)dx 
0 


1  h 


m2 


2h 
/         f(x)dx 

h 


Adding  terms  gives:! 


nh 
n  /    f(x)dx 

I      m.   = r- (2.11) 

i=l   x  n 


We  can  replace  in  (2.10)  to  obtain: 


x 

n 

Y    =    /    f(x)dx  +  Y  (2.12) 

0 


since  x  =  nh .   Equations  (2.10)  and  (2.12)  constitute  two  very- 
important  results.   Eq.  (2.12)  is  indeed  not  an  unexpected  one, 
and  several  considerations  can  be  drawn  from  it.   First  of  all 
the  value  of  Y   does  not  depend  on  the  step  size  h,  at  least 
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for  this  first  integration,  so  its  accuracy  is  determined  only 
by  the  exactness  of  the  integration  routine  used.   Second,  the 
values  of  Y(x)  can  be  directly  determine  at  any  arbitrary  point 
n,  that  is  to  say  we  need  not  find  its  previous  values  as  is 
the  case  for  many  numerical  methods.   Later  this  fact  will 
become  more  useful.   Finally,  we  see  that,  by  transposing  terms 
we  can  find  either  the  left  or  right  B.C.,  and  this  apparently 
unimportant  fact  is  indeed  a  key  step,  since  as  we  will  see 
later,  we  will  be  able  to  correlate  B.C. 's  of  higher  order 
differential  equations  and  solve  for  them  without  having  to 
generate  all  previous  values  of  the  unknown  function.   In  other 
words,  we  could  transform  a  Boundary  Value  Problem  into  an 
Initial  Value  Problem  and  vice  versa,  depending  on  what  is 
known,  and  what  we  are  looking  for. 


16 


III.   HIGHER  ORDER  DIFFERENTIAL  EQUATIONS 

A.   GENERAL 

The  previous  chapter  was  dealing  basically  with  the  main 
ideas  of  the  method,  and  we  have  solved  a  first  order  differ- 
ential equation.   We  are  now  going  to  extend  the  method  to 
higher  order  equations.   Essentially  what  we  will  do  is,  inte- 
grate several  times  the  right  hand  side  of  the  given  equation; 
as  required  for  its  degree,  using  the  same  approach  as  before. 
However,  we  should  keep  in  mind  that  it  is  only  for  the  first 
integration  that  we  will  use  the  given  function  f(x).   After 
this  integration,  this  function  is  no  longer  available  because 
what  we  have  is  a  set  of  straight  lines,  each  applying  to  a 
specific  interval.   We  need  to  find  now  a  method  that  will 
allow  us  to  perform  a  second  integration  of  the  function  f(x), 
from  the  given  set  of  lines. 

In  order  to  do  this,  let  us  consider  the  curves  of  Fig. 
3.1.  Suppose  we  have  been  given  the  following  second  order 
differential  equation: 


Y"(x)   =   f2(x)  (3.1) 


represented  in  Fig.  3 . 1  by  the  upper  curve.   Let  us  assume 


we 


have  determined  Y'(x)  =  f-^x)  after  a  first  integration  by  the 
method  just  introduced,  Let 
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re*) 


V»-f, 


Cx) 


i5?-^*+ 


y'cx)-  f, 


Cx) 


Vex)  =  f  ex) 


FIGURE     3.1 
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Y(X)   =   f(x) 

be  the  exact  result,  which  is  the  lower  curve  in  the  same 
figure.   Now  let 

y  '  (x)   =   m-,  x  +  c. 


and 


y ( x )   =   mx  +  c 

be  the  approximating  lines  to  f, (x)  and  f (x) ,  respectively, 

that  applies  inside  the  interval  from  a  to  b  only.   In  this 

figure,  we  can  easily  see  that  the  area  enclosed  by  the  upper 

straight  line  is  almost  equal  to  the  area  enclosed  by  the 

function  f, (x) .   Obviously  as  the  interval  becomes  smaller  and 

smaller,  both  areas  tend  to  be  equal;  in  the  limit  they  are 

indeed  equal.   It  is  evident  now  that  we  can  approximate  the 

exact  integration  of  f , (x) ,  by  integrating  the  line 

y'(x)  =  m, x  +  c, ,    provided  the  interval  is  small  enough  so  the 

error  is  negligible.   Now,  as  before,  we  need  to  determine  m 

and  c.   Recall  that  m,  and  c.  have  already  been  found.   So  by 

virtue  of  (2.5): 

b 
/   f1(x)dx 


m   = 


a 


h 
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but  by  (2.4) 


a 
m   =   — 


b 
/  (m,  x  +  c,  )  dx 


or: 


m 
m  =  -y(a  +  b)  +  c1  (3.2) 


So  we  know  now  the  slope  of  the  next  line;  the  line  belong- 
ing to  the  solution  function.   In  order  to  find  the  intercept 
c,  another  auxiliary  condition  must  be  supplied.   Let: 

Y(a)   =   Ya 

then,  replacing  this  condition  in  the  line  y(x)  =  mx  +  c,  we 
get : 


Y(a)   =   Y    =   ma  +  c 
a 


It  can  be  shown  that: 


y(x)   =   m(x-a)  +  Y  (3.3) 

a 


where  m  is  given  by  (3.2)  .   Eq .  (3.3)  gives  the  approximate 
solution  of  the  original  second  order  differential  equation, 
inside  the  interval  a  to  b.   Note  that  both  auxiliary  conditions 
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have  been  given  at  the  left  hand  side.   However,  it  really 
need  not  be  like  this,  as  will  be  shown  later. 

At  this  point  some  considerations  are  in  order.   With 
reference  to  Fig.  3.1  where  the  upper  curve  is  the  given 
function  Y"  (x)  =  f  2  (x)  ,  the  middle  curve  is  Y'(x)  =  f,(x), 
and  the  lower  one  is  the  solution  Y(x)  =  f(x). 

As  explained  previously  we  have  approximated  the  solution 
by  straight  lines.   Consider  now  the  slope  of  the  lower 
straight  line  m  which  is  given  by  (3.2)  .   This  equation  can  be 
written  as : 


m  =   -^[(irua  +  c,  )  +  (nub  +  c,  )  ] 


m   =   2[yl(a)  +  yl(b) ] 


m   =   -|[Ya  +  YbJ  (3.4) 


Since  the  ordinate  at  any  point  x  in  the  Y'(x)  curve  is 
indeed  the  slope  of  the  Y(x)  curve  at  the  same  point,  then 
by  virtue  of  (3.4) ,  we  can  say  that  the  slope  m  of  the  lower 
straight  line  is  the  average  of  the  slopes  corresponding  to 
its  end  points,  namely  point  Y   and  Y,  in  Fig.  3.1.   Recall 
that  these  two  slopes  are  given  exactly  by  the  ordinates  of 
the  middle  curve  Y'(x) .   Furthermore  let  us  evaluate 
y'(x)  =  m,  x  +  c-,  ,  the  equation  of  the  upper  straight  line,  at 
the  midpoint  of  the  interval,  that  is  we  let: 
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x   =  |(a  +  b) 


in  the  line  y' (x)  and  we  get: 


,  ,a  +  bv        ,a  +  b. 
y  ( — 2 — '      =  mi ' — 2 — '     1   =  m 


As  seen,  it  reproduces  m,  the  slope  of  the  lower  straight 
line . 

Now  we  can  have  a  deeper  insight  of  the  whole  process, 
and  why  it  works.   See  that  what  we  are  really  going  to  do  is 
to  use  the  slope  at  the  midpoint  (the  average  of  the  slopes  at 
the  extremes)  to  approximate  the  true  slope  of  the  straight 
line  y(x),  and  later  determine  the  corresponding  intercepts. 
In  fact,  had  we  known  exactly  the  function  Y'(x)  =  f, (x) ,  we 
would  have  been  able  to  determine  the  true  slope  of  y(x)  by 
direct  integration  as  before.   As  we  can  see  now,  we  are  intro- 
ducing a  source  of  error,  and  from  now  on  we  will  not  have  as 
accurate  results  as  in  the  first  integration. 

B.   DEVELOPMENT  OF  THE  METHOD 

Before  we  go  into  a  general  development  of  this  method  a 
slightly  different  notation  must  be  introduced,  and  we  will 
be  using  it  throughout  the  rest  of  this  research.   The  type 
of  differential  equations  we  will  be  looking  at  are  of  the 
form: 


Ylv(x)   =   f(x)  (3.5) 
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Such  equations  arise  in  the  study  of  deflections  of  beams, 
and  we  shall  apply  the  method  to  practical  beam  problems.   Let 
us  introduce  the  notation  to  be  followed.   We  will  write  the 
approximate  solution  line  of  (3.5)  as: 

y ( x )   =   mx  +  c 

in  a  clean  way  without  subscripts.   The  approximation  line 
to  the  first  derivative  will  then  be: 

y '  (x)   =   m,  x  +  c, 

and  for  the  second  derivative: 

y "  (x)   =   irux  +  c2 

and  so  on. 

The  boundary  conditions  will  be  represented  as: 


Y(xQ)   -   YQ         Y(xn)   =   Yn 


Y' (x  )   =   Y'        Y' (x  )   =   Y' 
o       o  n       n 


Y-(x  )   =   Y;        Y"(xn)   -   Y» 


etc.,  where  the  subscript  o  refers  to  the  left  hand  side,  and 
the  subscript  n  to  the  right,  of  the  range  under  consideration 
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There  is  another  subscript  that  we  must  keep  track  of. 
Recall  that  the  whole  solution  function  is  going  to  be  repre- 
sented by  a  set  of  approximating  lines,  see  Fig.  2.2.   Then  if 
we  start  at  the  left  end,  we  shall  use  the  subscript  n  as  a 
second  indice  to  identify  each  line  as  we  move  rightwards.   Be 
aware  that  we  will  take  the  subscript  i  immediately  to  the 
right  of  each  line  to  name  it  locally. 

Some  words  must  be  said  about  the  auxiliary  conditions. 
Recall  that  in  a  specific  problem  the  number  of  these  condi- 
tions are  equal  to  the  order  of  the  differential  equation  to 
be  solved.   If  all  these  conditions  are  given  at  the  same 
point,  then  we  are  dealing  with  an  Initial  Value  Problem. 
However,  if  the  auxiliary  conditions  are  given  at  both  ends 
of  the  range  of  interest,  then  we  have  a  Boundary  Value  Prob- 
lem.  When  we  are  able  to  have  an  exact  solution,  we  usually 
are  led  to  a  set  of  simultaneous  equations  involving  the  unknown 
B.C.'s;  those  can  be  solved  algebraically  and  lead  to  the 
whole  solution.   However  for  equations  with  no  exact  solution, 
like  the  ones  we  are  interested  in,  there  is  no  way  to  deter- 
mine the  conditions  at  the  other  end  of  the  interval,  and  that 
is  precisely  the  goal  of  our  approach.   We  need  to  correlate 
in  an  explicit  way  the  known  conditions  and  the  unknown  ones 
by  means  of  equations  of  the  form  of  (2.10) ,  so  we  can  be  able 

to  find  them. 

n 
Note  that  in  this  equation  the  term  h   £  m.  is  independent 

i=l  1 
of  B.C.'s  and  it  can  be  readily  found  even  if  we  do  not  know 
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them.   By  similarity  to  an  exact  problem  we  can  look  at  this 
term  as  being  the  general  solution  of  the  differential  equation. 
Note  that  this  term  actually  links  the  left  and  the  right 
boundary  conditions  in  a  very  explicit  way.   This  is  the 
heart  of  what  the  present  method  is  all  about.   With  this 
idea  in  mind  we  are  now  going  to  find  several  relationships 
for  successive  integrations. 

Let  us  start  by  noticing  that  (2.10)  is  a  general  relation 
that  applies  not  only  to  the  first  integration  but  also  to 
successive  ones.   The  only  term  that  is  particular  for  each 
integration  is  the  summation  of  the  slopes.   So  in  order  to 
perform  the  next  integration  we  shall  determine  this  term 
specifically  and  substitute  it  into  (2.10).   We  are  now  dealing 
with  a  fourth  order  differential  equation  since  this  type 
occurs  in  beam  problems.   The  application  to  other  orders  will 
be  self-evident.   Let 


y'V(x)   =  EA(x)  (3.6) 


After  a  first  integration  we  have: 


n 

Y"'   =   h   J   m.  +  Y,M  (3.7) 

n  ^3   l     o 

i=l 


where  the  subscript  3  in  the  summation  indicates  that  these 
slopes  belong  to  Y"*  =  f 3 (x)  and  are  given  by  (2.11).   We  need 
to  find  now 
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n 

Y"   =   h   J.   m.  +  Y"  (3.8) 

n         L2         1     o 

i=l 


Following  our  convention,  m0    is  given  by: 


m0     =   — ^(x.  ,  +  x.)  +  b-,  . 
2,i       2    l-l     i      3,i 


but 


3 ,i       l-l     3 , l  l-l 


then 


ra0  .   =   — 4 — (x.  ,  +  x.)  +  YV'   -  nu  •  x.  n 
2,i       2    l-l     l      l-l     3,i  l-l 


or 


m2,i   "   ^h  +  Xi-1  (3-9) 


Now,  since 


n-1 

:"'  ,   =   h   I.  m. 
n-1        L3      l 

i=l 


then 
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m^  n-1 

m2/n   =   -^h  +  h   ^  m.  +  Y-  (3.10) 

i=l 


This  last  relationship  in  the  way  it  has  been  derived 
becomes  a  general  one,  and  expresses  the  slope  m-    at  any 
point  i  as  a  function  of  the  previous  slopes.   Similar  rela- 
tions can  be  derived  for  successive  integrations  just  by 
shifting  the  first  indice.   Now  we  need  to  find  the  summation 
of  these  slopes.   It  becomes: 


n  ,   n  n-1    i 

Y0  m.   =   £   Y,  m.  +  h   V     V   m.  +  nY"'  (3.11) 

i=l  i=l         1_1   j=l 


Now,    replacing    (3.11)    in    (3.8) 


Y»      =      h[£      Y,    m.    +   h      y         Y      m.    +   nY"'  ]    +   Y" 
n  2      ^3      i  ■    i         3      1  o  o 

i=l 


i=l  j=l 


or 


:-      =      h2[±      Lm,    +      I  I-   m,]    +   nhY"'     +   Y"  (3.12) 

n  2      ^3      l         .  L  -.         L3      i  o  o 


i=l 

1=1  :=i 


Let : 


n-1      i 


i=l  i=l      j=l 
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Note  that,  as  stated  above,  this  last  term  S2IS  also 
independent  of  boundary  conditions  and  can  be  readily  deter- 
mined.  Equation  (3.12)  gives  the  particular  solution  to  a 
second  order  differential  equation  as  well  as,  in  this 
specific  case,  it  represents  the  second  derivative  of  Y(x) . 

The  whole  process,  so  far,  can  be  summarized  as  follows: 
Equation  (2.10)  gives  the  solution  at  any  point  n  of  any 
integration,  first  second,  etc.   In  order  to  use  this  equation 
we  need  to  find  the  summation  term,  and  this  is  given  by  a 
relation  of  the  form  of  (3.11),  which  is  also  a  general  rela- 
tionship and  expresses  the  summation  of  slopes  in  terms  of 
the  previous  one.   Recall  that  the  only  slope  summation  we  know 
is  the  one  involving  the  known  function  f (x) ,  and  is  given  by 
(2.11).   That  is  why  all  further  summations  must  be  expressed 
as  a  function  of  this  one.   This  last  summation  is  then  re- 
placed in  (3.8),  and  simplified  if  possible.   The  next  step  is 
to  identify  the  terms  that  are  independent  of  boundary  condi- 
tions and  isolate  them  as  in  (3.12)  and  to  evalute  them 
separately. 

We  can  continue  in  the  same  way  until  we  get  a  solution 
formula  for  Y(x).   It  can  be  shown  that  the  complete  solution 
for  a  fourth  order  differential  equation  of  the  form  Y  (x)  =  f (x) 
is  given  by  the  following  relationships: 


Y'"   =   S-,  +  Y"1  (3.13) 

n        3     o 


Y"   =   hSn  +  nhY,M  +  Y"  (3.14) 

n        2       o      o 
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Y'       =      h2S,     +    N1h2Y"'     +    nhY"    + 
n  1  1         o  o 


,     +    N,h    Y"'     +    nhY"    +   Y'  (3    15) 

1  1        o  o  o  v  ; 


Yn      =      h3S^    +    N   h3Y,M     +   N,h2Y"    +    nhY'    +    Y  (3.16) 

n  ooo  lo  oo  ' 


where : 


n 
S3      =         I      m.  (3.17) 

i  =  l      1 


-■       n  n-1      i 

S2      =      j      I      m±   +      I         I      m  (3.18) 

i=l  i=l    j=l      J 


n  n-1      i  n-2      i 


l  2  n~z      -1        2 

!1      =      J      I      mi   +      I         I      mi    +      I         I         I      mk  (3.19) 

i=l  i=l    j=l      J         i=l    j=l    k=l 


-,       n  ~    n-1      i  -.    n-2      i         j  n-3      i         j         k 

;      =    o       J   m-    +   t      J         J    m-    +   4      7         I         L    m,     +      I         1         1         I   mn 
i=l  i=l    j=l    J  i  =  l    j=l    k=l  i=l    3  =  1    k=l   1=1 

(3.20) 


and 


n-1 
N,       =      |+      I         i  (3.21) 

1  i=l 


n-1  n-2      i 

N         =      j  +       I  i    +       I  I       j  (3.22) 

i=l  i=l      j 
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m.       = 

1 


1 

h 


(i-Dh 


ih 
/        f(x)dx 


1    <    n 
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IV.   ERROR  ANALYSIS  AND  CONVERGENCE 

Refer  to  Fig.  4.1  where  we  have  plotted  the  functions 
Y1  =  f,  (x)  and  Y  =  f (x)  together  with  their  respective 
approximating  straight  lines  inside  the  interval  x  =  a  to 
x  =  b.   We  assume  that  the  upper  straight  line  belongs  to  a 
first  integration;  in  other  words  we  have  not  introduced  any 
error  so  far.   A  second  integration,  however,  will  be  performed 
using  this  line  y, (x)  since  we  do  not  know  the  function  f , (x) . 
It  is  worth  considering  now  what  happens  when  we  integrate 
the  line  instead  of  the  function  itself.   Let  m  be  the  slope 
of  the  lower  straight  line.   This  is  given  by: 

b 

/   f  (x)dx 

a 
m  =   — 


b-a 


Let  m*  be  the  approximating  slope  given  by: 

b 
j       (m,  x  +  c,  )  dx 


m*   =   —  (4.1) 

b-a 


Now  for  the  particular  case  shown  in  the  figure,  we  have: 


b  b 

/   (mxx  +  c1)dx   >    /   f1(x)dx         (4.2) 
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FIGURE    4.1 


32 


since  the  area  enclosed  by  the  curve  is  smaller  than  the  area 
enclosed  by  the  line  (we  assume  that  there  are  no  inflection 
points  and/or  discontinuities  inside  the  interval) .   It  is 
now  clear  that: 

m*   >   m 
and  the  error  introduced  is: 


b  b 

e    =    /  (m1x+c;.)dx  -   /   f, (x)dx      (4.3) 


a-b 


which  is  equal  to  the  shaded  area  in  Fig.  4.1. 

As  indicated  in  Chapter  III,  m*  comes  out  to  be  the  average 
of  the  slopes  at  the  end  points  of  the  lower  straight  line, 
and  it  is  this  slope  we  work  with.   So,  the  actual  approximat- 
ing line  y*(x)  =  m*x  +  c*,  shown  in  Fig.  4.2  is  steeper  than 
y(x)  since  m*  >  m.   The  distance  ED  becomes  the  error  introduced 
and  we  can  not  evaluate  it  because  we  do  not  know  f (x) .   Another 
integration  using  y* (x)  will  make  things  even  worse.   Note 
that  the  error  to  be  committed  this  time  will  be  larger  and 
is  given  by  the  shaded  area  in  Fig.  4.2. 

So  far  it  appears  that  the  method  is  divergent  in  nature 
due  to  the  fact  that  the  error  will  grow  bigger  and  bigger 
unless  we  find  some  means  to  correct  it.   At  this  point,  the 
only  possible  way  to  do  this  is  by  reducing  the  "shaded" 
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areas,  which  means  that  we  have  to  reduce  the  step  size. 
Fortunately,  we  can  say,  at  least  in  theory,  that  in  the  limit 
there  will  be  no  error;  in  practice  however,  a  very  small 
step  size  will  increase  the  round-off  errors.   The  examples 
provided  show  that  the  approximation  is  indeed  acceptable. 

There  is  a  special  case,  however,  where  an  error  correc- 
tion can  be  made.   This  applies  only  to  boundary  value  problems 
Recall  that  a  solution  of  any  differential  equation  by  this 
method  is  given  by  the  general  relationship: 


n 

Y    =   h   J   m.  +  Y 
n        .  L  -.       1    o 
i=l 


In  B.V.P.,  we  know  the  conditions  at  both  ends  of  the  interval 
So  this  last  equation  can  be  written  as: 


n         Y  -  Y 

I      m.  =   -^-^  (4.4) 

.  L ,       l        h 
i=l 


where  the  right  hand  side  is  known.  Now  since  we  are  actually 
dealing  with  approximations  to  the  true  slopes,  what  we  really 
have  as  a  summation  is  the  next  term,  call  it  M: 


n 

M   =    y   m* 
i=l   X 
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Our  approximate  solution  is  then: 


n 
Y*   =   h  I      m*  +  Yq  (4.5) 

i=l 


and  the  error  at  point  n  will  be  given  by 


e    =   Y*  -  Y  (4.6) 

n       n     n 


Now  let  us  see  how  the  correction  could  be  carried  out 
individually  at  each  interval.   To  do  this  we  must  find  the 
exact  slope  m.   It  is  given  by: 

m  =   m*  -  e 

where  e  is  the  shaded  area  of  Fig.  4.1  and  is  given  by  (4.3) 
The  slope  summation  is  therefore: 


*  -   T 


)   m.   =    >   m?  -  ne 
i=l   x      i=l   => 


Multiplying  by  h  and  adding  Y   to  both  sides,  we  obtain 


htm.  +  Y    =   h   )   m*  +  Y   -  ihe 
.  L,   i     o        .  L1       1     o    J 
i=l  i=l 


or 
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Yj    =    Yj  "  jhe  (4-7) 


When  j  =  n,  we  have 


Y    =   y*  -  nhe 
n       n 


or 


Y*  -  Y    =   nhe 
n     n 


but  from  (4.6)  we  get 


e  =  HE  <4-8' 


Replacing  this  last  relation  in  (4.7) 


Y  .   =   Y*  -  1  e 
:       3    n   n 


and  the  general  relationship  given  by  (2.10)  transforms  into 


j 

Y.   =hTm*+Y-ie  (4.9) 

j        .  L ,       i  o    n   n 

J        i=l 


which  is  the  corrected  solution  at  any  point  j.  In  summary, 
what  we  need  to  do  to  correct  the  solution  is:  a)  determine 
Y*  by  performing  the  integration  without  any  correction  by 
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using  equation  (4.5);  b)  compare  Y*  with  the  given  B.C.  Y 

and  obtain  e;  c)  Perform  again  the  integration  this  time  using 

the  relationship  given  by  (4.9). 

From  (4.9)  we  can  find  a  whole  set  of  corrected  equations 
similar  to  those  given  by  (3.13)  to  (3.20).   However,  by  doing 
this  we  will  enormously  complicate  those  relationships  and 
the  computational  effort,  together  with  round-off  errors,  may 
not  give  any  advantage  at  all,  especially  since  we  expect  the 
local  error  e  to  be  very  small. 

There  is,  however,  another  stronger  reason  not  to  do  that. 
The  fact  is  that  we  have  assumed  that  the  error  e  is  a  con- 
stant and  this  is  not  always  the  case.   If  we  restrict  our- 
selves to  integrate  linear  and/or  constant  functions,  then  the 
corrected  method  could  be  justified  since  for  these  functions 
the  error  given  by  (4.3)  is  a  constant,  for  the  shaded  area 
of  Fig.  4.1  will  always  be  the  same.   In  more  general  cases, 
the  error  e  will  be  unpredictable  and  no  correction  can  be 
made . 

This  example  serves  only  to  illustrate  how  the  error 
correction  could  be  carried  out.   But  it  would  not  be  practical 
and  a  reduction  in  step  size  appears  to  be  the  best  correction 
as  we  will  see  in  the  examples. 

As  stated  above,  the  method  appears  to  be  divergent  in 
nature,  so  our  approximate  solution  will  look  like  Fig.  4.3. 
In  this  figure,  note  that  for  any  integration  after  the  first 
one,  the  set  of  approximating  lines  diverge  from  the  exact 
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solution,  the  further  we  move  to  the  right,  the  greater  the 
error  introduced.   In  dealing  with  I.V.P.'s,  the  points  Y1 
and  Y   are  known  and  these  are  our  starting  conditions.   The 
solution  will  appear  then  as  shown  in  this  figure.   Since  we 
do  not  know  the  exact  end  points  at  the  right  hand  side, 
namely  points  Y'  and  Y,  there  remains  an  uncertainty  in  the 
accuracy  of  the  solution. 

For  B.V.P.'s,  however,  the  other  extreme  points  are  known. 
Referring  to  Fig.  4.3,  suppose  we  are  given  as  boundary  condi- 
tions, points  Y'  and  Y   in  the  lower  curve  of  this  figure. 
The  method  we  are  dealing  with  requires  that  we  know  the 
initial  points  to  be  able  to  start  the  integrating  algorithm. 
If  we  know  the  exact  starting  points,  namely  Y'  and  Y  ,  then 
by  using  these  initial  conditions,  our  solution  will  appear 
as  shown  in  Fig.  4.3.   That  is  to  say  we  will  not  end  up  at 
point  Y  which  is  exact,  but  at  point  Y*.   To  be  able  to  reach 
the  exact  point  Y,  which  we  assume  is  known,  we  need  to  give 

the  algorithm  a  "wrong"  starting  point  Y'*.   It  is  this  approxi- 

o 

mate  starting  point  which  we  find  using  the  relationships  given 
by  equations  (3.13)  to  (3.22)  that  allow  us  to  match  the  exact 
solution  in  the  whole  range  as  shown  in  Fig.  4.4  in  the  lower 
curve . 

Note  that  if  we  were  using  corrected  relationships  we 
would  be  able  to  supply  the  algorithm  with  "exact"  starting 
values,  and  still  match  the  exact  end  points.   Unfortunately, 
it  is  now  the  upper  set  of  solution  lines  that  absorb  the 
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inherent  error  of  the  method  since  we  have  given  them  approxi- 
mate initial  values.   The  solution  now  is  going  to  look  like 
Fig.  4.4.   The  lower  solution  becomes  "exact"  but  the  upper 
actually  shifts  apart  from  the  exact.   We  can  explain  this  by 
saying  that  the  known  B.C.'s  are  fixed  by  the  algorithm  while 
the  free  (unknown)  ones  will  take  the  error,  and  this  is 
exactly  what  happens  in  the  examples  in  the  next  section.   In 
summary,  all  we  do  is  to  shift  the  error  from  one  place  to 
another.   We  can  not  get  rid  of  it  unless  we  use  corrected 
relationships  or  an  infinite  number  of  intervals  to  minimize 
the  local  error. 
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V.   APPLICATIONS  AND  RESULTS 

A.   GENERAL 

In  this  chapter  we  are  going  to  study  a  particular  appli- 
cation of  the  integrating  algorithm  introduced  in  the  previous 
sections,  namely,  the  deflection  of  beams,  where  a  fourth 
order  linear  differential  equation  occurs,  and  several  combina- 
tions of  B.C.'s  can  be  considered. 

From  Mechanics  of  Solids,  the  relationships  governing  the 
deflection  of  beams  are  given  by: 


i  v 
EIY   (x)   =   q(x)      (load  intensity)        (5.1) 


EIY"'  (x)   =   V(x)      (shear  force)  (5.2) 

EIY" (x)   =   M(x)       (moment)  (5.3) 

Y ' (x)   =   slope 

y(x)   =   deflection 

The  method  we  are  dealing  with  is  basically  a  method  of 
successive  integration.   Because  of  that  we  can  start  to  inte- 
grate from  any  of  the  above  relationships  provided  we  know 
explicitly  the  right  hand  side  of  the  equation.   In  some  cases, 


43 


a  given  problem  could  be  solved  by  integrating  the  moment 
second  order  differential  equation.   By  doing  this,  we  would 
be  able  to  achieve  more  accuracy  since  only  two  integrations 
need  to  be  performed  and  the  accumulated  error  will  be  rela- 
tively small. 

However,  when  the  loading  of  the  beam  is  a  complicated 
distribution  and  the  expression  for  the  bending  moment  is 
difficult  to  obtain,  then  the  fourth  order  differential  equa- 
tion will  be  the  one  to  use.   One  advantage  in  using  this 
expression  is  that  the  integrating  algorithm  will  provide 
results  for  the  shear  force,  bending  moment,  slope  and 
deflection  simultaneously.   When  using  the  second  order  equa- 
tion we  only  get  slope  and  deflection. 

Since  eq.  (5.1)  is  a  more  general  relationship,  we  are 
going  to  use  this  expression  in  the  examples.   In  fact,  most 
problems  can  be  expressed  in  that  way.   Recall  that  a  concen- 
trated force  can  be  treated  as  a  distributed  force  acting  in 
a  very  small  interval,  and  a  concentrated  moment  reduces  to  a 
couple  of  concentrated  loads  acting  in  opposite  directions  a 
small  distance  apart.   Example  2  illustrates  this  procedure. 

The  examples  provided  in  this  chapter  account  for  all 
types  of  B.C.'s.   Fig.  5.1  shows  the  4  cases  to  be  treated, 
and  we  shall  study  them  one  by  one  in  the  next  subsections. 
In  all  cases  the  flexural  rigidity  EI  and  the  length  L  is 
set  equal  to  1,  for  simplicity.   Discontinuous  types  of  loading 
are  emphasized,  and  of  course,  many  combinations  of  loadings 
can  be  dealt  with  by  simply  using  the  principle  of  superposition. 
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The  computer  program  used  to  solve  these  problems  is  shown 
at  the  end  of  this  chapter,  and  it  is  a  straightforward  coding 
of  equations  (3.13)  to  (3.22) .   This  programming  technique,  in 
particular,  tries  to  minimize  round-off  errors.   For  now,  it 
suffices  to  say  that  all  problems  are  solved  using  SINGLE 
PRECISION,  one  hundred  intervals  (problem  2  uses  200),  and  we 
get  accuracy  to  the  fourth  and  even  to  the  fifth  decimal  places 
in  some  examples.   Finally,  only  the  fourth  integration,  the 
deflection,  is  compared  with  the  exact  solution.   Recall  that 
the  error  is  cumulative,  and  it  is  here  that  we  expect  the 
bigger  error.   We  may  conclude  then  that  preceding  integrations 
are  more  exact.   However,  we  should  keep  in  mind  that  an 
approximated  starting  point  is  needed  in  some  cases,  and  the 
results  will  be  shifted  by  some  amount  from  the  exact.   We 
will  see  this  as  we  proceed  into  the  next  section. 

B .   EXAMPLES 

1.  Example  1;   Clamped-Free  Beam 

Refer  to  Fig.  5.1a.   This  is  a  statically  determinate 
cantilever  beam  loaded  over  half  its  length  by  a  uniformly 
distributed  load  of  1  unit  of  weight  per  unit  of  length  as 
shown.   We  shall  determine  the  shearing  force,  bending  moment, 
slope  and  deflection  at  discrete  points  of  the  beam. 

2 .  Solution 

This  is  a  B.V.P.  but  for  purpose  of  illustrating  how 
to  use  the  method  when  dealing  with  I.V.P.'s,  we  will  treat 
it  as  such.   All  the  required  initial  conditions  can  be  drawn 
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from  the  given  configuration.   Here  the  function  f (x)  to  be 
integrated  is  discontinuous  and  given  by  a  set  of  two  differ- 
ent functions  defined  by: 


f1(x)   =   0 


0  <  x  <  0.5 


f2(x)   =   1.0 


0.5  <  x  <  1.0 


We  need  to  determine  the  initial  values  to  start  the 
algorithm.   From  statics,  it  can  easily  be  shown  that  the 
initial  conditions  are  as  follows: 


Shear : 

Y  "• 

o 

= 

-0.5 

n 

— 

0.0 

Moment: 

Y" 
o 

= 

0.375 

Y" 
n 

= 

0.0 

Slope : 

Y' 

o 

= 

0.0 

Y' 
n 

= 

Unknown 

Deflect 

ion: 

Y 
o 

= 

0.0 

Y 
n 

= 

Unknown 

By  supplying  these  initial  values  to  the  program,  we 
obtain  the  computer  output  shown  on  the  next  page.   The  last 
column  shows  the  exact  deflections.   In  these  results,  note 
that  all  the  initial  values  are  exact,  and  that  the  last  value 
found,  namely  the  deflection  of  the  beam  at  x  =  1.0,  is  expected 
to  have  the  biggest  cumulative  error.   However,  as  we  can 
see,  we  get  3  to  4  digits  accuracy.   Here  we  can  conclude 
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EXAMPLE    1      . 
CLAMPliC-FREE    EcMA 


BCUfsGAF-Y 

CONCITICNS: 

Y(0)  = 

O.CCOO 

Y     (C)=       O.JOOO 

tf(0)  = 

0.3750 

V(0  )=    -0.5000 

X: 

SHE4R 

^cmsnt 

SL3P3 

DEFLECTION 

EXACT    DEF 

0.00 

-o.scccoo 

0.375000 

O.COOOOO 

C. 000000 

-0.000000 

0.05 

-O.50C0OO 

G. 3 5 000  0 

0  .C18125 

C. 000453 

C. 000458 

0.10 

-o.5:ccoo 

C. 325000 

U.C35000 

C.0j1792 

0.001792 

0.15 

-o.5accoo 

C. 300000 

O.C 50625 

C. 003937 

0.003933 

0.2G 

-0.500000 

Q. 275000 

0  .C65000 

C.  006333 

0.006333 

0.25 

-C.5CCCG0 

G.25GG0G 

0.C73125 

C. 010416 

0.010417 

0.3C 

-C.5CCC0O 

0.225CJ0 

O.C90000 

C. 014625 

0.014625 

0.35 

-0.  5CCCG0 

0.200000 

0.  100625 

C. 01 9395 

0.019396 

0.40 

-0.5CCC00 

0.175000 

O.llOuOO 

C.02466o 

0.024667 

0.45 

-o. 5:ccoo 

C.15uOOG 

0.116125 

C. 0^0375 

0.030375 

0.50 

-o.5:cooo 

C. 125000 

0.12D0OO 

C. 036453 

0.036459 

0.55 

-0.43CCJ0 

0.10125  0 

0  .130646 

C.04285h- 

0.042355 

O.oG 

-O.43C0O0 

G.C80000 

0.135167 

C. 049504 

0.C49505 

0.65 

—  0.3  5CC00 

0.061250 

0  .138638 

C. 05b35t 

0. 056355 

0.7C 

-0.3CCCO0 

Q.C+50Q0 

0  .141334 

C. 063353 

0.063358 

0.  75 

-G.25CC00 

O.OJ1250 

0.143230 

C. 070475 

0.070475 

O.dC 

-0.2CCC00 

C.02G00O 

0.144501 

C.J77o70 

0.077o71 

0.35 

-0.150000 

C.  011250 

0  .145272 

C. 034916 

0.034916 

0.90 

-o.iaccoo 

C.CC5000 

0.1456o7 

C. 092191 

0.092191 

0.95 

-0. C5CCG0 

(J.CG1249 

0.145313 

C. 099479 

0.099479 

l.GC 

-o.oacooo 

-C.U00001 

0.  145334 

C. 106770 

0.106770 
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that  previous  values  of  all  other  columns  are  even  more 
accurate . 

3.  Example  2:   Clamped-Roller  Beam 

Refer  to  Fig.  5.1b.   This  is  a  statically  indeterminate 
beam  loaded  by  a  moment  of  10  units  of  weight  x  unit  of  length, 
acting  at  the  middle  of  the  beam.   As  before  we  are  going  to 
determine  the  shear,  moment,  slope  and  deflection  of  the  beam. 

4 .  Solution 

This  example  will  emphasize  three  things:   How  to  deal 
with  a)  concentrated  forces,  b)  concentrated  moments  and  c)  the 
necessary  initial  conditions  required  for  this  particular 
case . 

a.   Concentrated  Forces 

We  can  conceive  of  a  concentrated  force  as  a  dis- 
tributed force  of  high  intensity  distributed  over  a  very  small 
length  of  beam.   Recall  from  the  previous  chapters  that  the 
slope  of  any  approximating  line  at  any  point  i  is  given  by 
(2.9)  where  f(x)  is  the  distributed  load  for  the  case  of  beams. 
In  this  equation  the  integral  in  the  numerator  represents  the 
area  under  the  f(x)  curve.   This  area  is  actually  the  total 
load  in  the  specified  interval.   By  decreasing  the  interval 
while  still  holding  the  same  inside  area  constant  (same  total 
load) ,  the  distributed  force  tends  to  a  concentrated  force  of 
equal  magnitude  as  the  total  load.   In  the  limit  it  becomes  a 
concentrated  force  acting  on  a  mathematical  point  on  the  beam. 
For  our  purpose,  in  the  integrating  algorithm  we  simply  let: 
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m.   =   0      for     x  4   x. 
1  i 


P      r 
m.   =  r-     for     x  =  x. 
l     h  i 


/  f  (x)  dx   =   P      for     x  =  x. 

where  P  is  the  concentrated  force  acting  at  point  x.. 
b.   Concentrated  Moments 

A  concentrated  moment  can  be  treated  in  a  similar 
way  if  we  use  the  second  order  differential  equation.   However, 
in  our  case  we  need  to  decompose  the  given  moment  into  a  pair 
of  concentrated  forces  of  equal  magnitude  acting  in  opposite 
directions  at  equidistant  points  from  the  location  of  the 
given  moment.   From  statics  we  know  that: 

Moment   =   Force  x  Distance 

Here  we  can  have  several  combinations  of  force 
times  distance  provided  we  keep  this  product  a  constant  equal 
to  the  given  concentrated  moment.   However  as  we  will  see 
later,  accuracy  is  achieved  only  when  we  use  very  small  dis- 
tances and  consequently  very  large  concentrated  forces.   In 
this  way  of  moment  decomposition,  the  usual  behavior  of  a 
concentrated  moment  is  achieved.   For  larger  distances  the 
couple  is  not  an  accurate  representation  of  a  concentrated 
moment  and  the  results  are  not  so  accurate. 
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In  the  current  example  we  use  M  =  10  and  it  is 
decomposed  as  follows: 

10   =   1000  x  0.01 

a  pair  of  opposite  forces  of  1000  units  acting  0.01  units  of 
length  apart  from  each  other.   In  our  problem  we  use  200 
intervals  of  0.00  5  each,  one  concentrated  force  is  located 
at  x  =  0.495,  and  the  other  opposite  force  of  equal  magnitude 
is  placed  at  x  =  0.50  5.   So  the  loading  function  to  be  inte- 
grated is  defined  by  5  partial  functions  as  follows: 


f1(x)   =   0 


0  <  x  <  0.495 


Jf9(x)dx   =   -1000 


x  =  0.495 


f3(x)   =   0 


0.495  <  x  <  0.505 


Jf.(x)dx   =   1000 


x  =  0.505 


f5(x)   =   0 


0.505  <  x  <  1.0 


c.   Necessary  Initial  Conditions 

We  need  to  consider  now  the  B.C.'s.   For  the  present 
configuration  we  have: 
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Shear : 

Y  "• 
o 

= 

Unknown 

Moment : 

Y" 
o 

= 

Unknown 

Slope : 

Y' 
o 

= 

0 

Deflect. 

ion: 

Y 
o 

= 

0 

Y  "  • 

n 

— 

Unknown 

Yn 

= 

0 

Y' 
n 

= 

Unknown 

Vn 

= 

0 

To  start  the  algorithm  we  need  to  determine  the  initial 

shearing  force  Y"1  and  the  initial  bending  moment  Y".   From 

o  o 

equations  (3.13)  to  (3.16),  after  replacing  the  known  values 
we  get: 


0   =   hSn  +  nhY"1  +  Y" 
2       o      o 


Y'   =   h2S,  +  N,h2Y"'  +  nhY" 
n        1     1    o       o 


0   =   h3S   +  N  h3Y,M  +  N,h2Y" 
o    o    o     1    o 


This  is  a  set  of  three  simultaneous  algebraic  equations  with 
three  unknowns  and  we  solve  for  the  initial  conditions.  The 
results  are: 


•  ii   


h(N  S0  -  nS  ) 
o  2 o_ 

nN.  -  N 
1    o 


(5.4) 


S   -  N,S~ 
v»i  =    °     1  2 
o      nN,  -  N 


(5.5) 


52 


All  the  terms  involved  in  these  expressions  are 
defined  by  the  relations  (3.13)  to  (3.22),  and  they  can  be 
computed  in  advance  since  they  are  independent  of  boundary 
conditions.   These  initial  values  are  only  approximations  as 
we  know,  but  serve  to  match  the  given  B.C.'s  at  the  other 
extreme.   The  results  are: 


Shear:  Y"'   =   11.21219  Exact   =   11.25 

o 


Moment:         Y"    =   -1.21219  Exact   =   -1.25 

o 


Running  the  program  with  these  initial  values,  we 
get  the  results  shown  on  the  next  page.   The  deflection  appears 
accurate  to  about  2  percent  for  the  more  meaningful  values  of 
deflection.   However,  if  we  reduce  the  distance  between  the 
opposite  forces,  we  improve  the  accuracy.   But  if  we  increase 
it,  the  discrepancies  grow  enormously  even  for  a  very  small 
increment,  and  it  does  not  reflect  a  concentrated  moment 
behavior. 

The  slope  shows  the  expected  behavior.   However,  we 
should  not  expect  high  accuracy  since  we  have  started  the 
algorithm  with  approximate  starting  points  of  shear  and  moment. 

At  x  =  1.0,  we  have  Y   =  -0.631.   The  exact  is  -0.625.   Pre- 

'  n 

ceding  values  of  slope  are  expected  to  be  better  since  here 

we  start  with  an  exact  known  initial  point.   Similarly  for 

the  moment  we  start  with  an  approximation  but  get  more  accurate 
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EXAMPLE    2      ' 
CLAMPtC-SGLLE*     EEAM 
BOUNDARY    CONCITICNS 


a:  SHt^R 

0.00  11.212190 

0.05  11.212190 

0.10  11.212190 

0.15  11.212190 

0.2C  11.212190 

0.25  11.212190 

0.3C  11.212190 

0.35  11.212190 

Q.4C  11.212190 

0.45  LI. 212190 

G.5U  -983.737500 

0.55  11.212190 

0.6C  11.212190 

0.6  5  11.212190 

0.7C  11.212190 

0.75  11.212190 

0.3C  il. 212190 

0.85  11.212190 

0.9C  11.212190 

0.95  11.212190 

1.00  11.212190 


YIO) =0.0 
Y    (0)=0.Q 


Y(  1J  =  0.C 
Ml  l)=J.O 


fCMENT 

1.212193 

■C. 651534 

•C.C9C975 

C.46S634 

1.020243 

1.59C853 

2.151463 

2.712071 

3.Z 72630 

3.833291 

■2.106097 

■5.045436 

-4.48<+873 

-2.S24269 

-2.363659 

-2.8CJ050 

-2.2424*0 

-1.6ol330 

-1.121222 

-C.560ol3 

-0.000003 


SLOPE 

C.COOOOO 

■0  .C46594 

■0.C65158 

■0.C55692 

-0.C13195 

0.C47333 

0  .140890 

0.262479 

0.412098 

0 .589747 

0.764177 

0.504137 

J .265877 

0.C55649 

-0.  126550 

-0.23U717 

-0.406855 

-0.504961 

-0.5750^7 

-0.617034 

-0.631098 


DcFLcCTION         EXACT    DEF 


C. 000000 

-C.0J123O 

-C. 004190 

-C.0J7327 

-C. 009290 

-C. 008677 

-C. 004087 

C. 005382 

C. 022630 

C. J4  756  1 

C. 081965 

C. 113693 

C. 132333 

C. 1*0755 

C. 138867 

C.j.23570 

C. 111265 

C. 033354 

C. 061239 

C. 031319 

C. 000000 


0.000000 

■0.301327 

-0.00^374 

-0.007734 

-0.009999 

-O.0U9766 

-0.005624 

O.00J823 

C.  020000 

G.o4<+297 

0.073125 

0.110390 

C.130COG 

0.138359 

0.13o875 

0.12o953 

0.110000 

0.03  7  422 

0.060625 

0.031016 

C.OOOCOO 
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as  we  proceed,  since  we  end  up  with  the  expected  result  of 
0.0  at  x  =  1.0 . 

In  the  case  of  the  shear  force,  we  see  it  is  con- 
stant throughout  the  beam  length  as  expected  but  it  is  shifted 
by  an  amount  of  0.0378  from  the  exact.   It  is  interesting  to 
note  that  at  the  point  of  the  loading  we  have  a  very  strong 
shearing  force.   For  a  pure  moment  acting  at  this  place,  we 
shall  neglect  this  result.   However,  if  we  do  have  an  actual 
couple  than  a  big  shearing  at  this  point  should  be  expected. 
Note  that  the  algebraic  sum  of  shear  forces  gives: 

11.21219  -  (-988.7875)   =   1000.0 

which  is  indeed  the  couple  acting  at  this  point. 

5.  Example  3:   Pin-Roller  Beam 

Refer  to  Fig.  5.1c.   This  is  a  statically  determinate 
simple  supported  beam  subjected  to  a  triangular  type  of  load 
distributed  over  the  central  portion  of  the  beam  as  indicated 
in  the  figure.   We  need  to  determine  the  shear,  moment,  slope 
and  deflection  of  the  beam. 

6 .  Solution 

First  of  all,  some  comments  about  this  configuration 
are  in  order.  This  example  was  chosen  in  order  to  emphasize 
the  ease  with  which  we  can  switch  from  one  kind  of  load  dis- 
tribution to  another.  As  explained  at  the  beginning  of  this 
work,  it  is  the  fact  that  the  algorithm  solves  a  "new"  differ- 
ential equation  at  every  single  step  that  allows  us  to  deal 


55 


with  different  loads.   Therefore  it  is  possible  to  have  a 
different  kind  of  loading  per  subdivision.   Another  reason  for 
selecting  this  problem  is  that  here  we  deal  with  linear  forcing 
functions  and  we  are  going  to  perform  four  successive  integra- 
tions of  a  linear  function.   This  will  lead  to  a  polynomial 
of  fifth  degree  and  a  bigger  cumulative  error  is  to  be  expected 

Proceding  with  the  solution,  the  distributed  load 
f(x)  to  be  integrated  is  given  by  four  different  functions: 


f1(x) 


=   0 


0  <  x  <  0.25 


f2(x) 


=   4x-l 


0.25  <  x  <  0.50 


f3(x) 


=   -4x+3 


0.50  <  x  <  0.75 


f4(x) 


=   0 


0.75  <  x  <  1.0 


As  indicated  above,  the  computer  program  is  such  that 
it  can  switch  from  one  kind  of  loading  to  the  other  at  the 
specified  point.   Now,  in  order  to  start  the  algorithm  we 
need  to  determine  the  initial  points  from  the  B.C.'s  of  a 
simple  supported  beam.   These  are: 


Shear: 

Y  "• 

o 

Unknown 

Y... 
n 

Unknown 

Moment : 

Y" 
o 

= 

0 

Y" 

n 

= 

0 

Slope : 

Y' 
o 

= 

Unknown 

Y' 

n 

= 

Unknown 

Deflect 

ion : 

Y 
o 

= 

0 

Y 
n 

= 

0 
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Now,  from  the  known  relationships  after  we  substitute  in  the 
given  values,  we  get: 


0   =   hS~  +  nhY,M 
2       o 


0   =   h3S   +  N  h3Y,M  + 


nhY' 
o    o    o      o 


This  is  a  pair  of  simultaneous  algebraic  equations  in  two 

unknowns.   We  solve  for  Y1"  and  Y1.   The  results  are: 

o       o 


Y'   =   (£)2(N  S9  -  nS  , 
o      n    o  2.  o 


(5.6) 


Y  "• 
O 


n 


(5.7) 


Similarly,  the  computer  supplies  the  values  of  the 
variables  involved  in  these  relations  and  these  are  as  follows 


Shear: 


O 


=   -0.12499 


Exact   =   -0.125 


Moment : 


Y' 
o 


=   0.014972 


Exact   =   0.014974 


which  shows  a  remarkable  accuracy.   Running  the  program  with 
these  initial  values  we  obtain  the  results  shown  on  the  next 
page.   A  comparison  between  the  approximate  deflection  and  the 
exact  reveals  an  error  of  about  2  percent  at  the  point  where 
the  largest  deflection  takes  place.   In  this  example,  we 
expected  a  larger  error,  but  the  algorithm  appears  to  perform 
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EXANPLE 

3     • 

PIN-ROLLER    3E^ 

BOUNCARY 

CONCniCNS  : 

Y(0) =0. 

.0            Yll)=0 

.0 

M(0)=0, 

.  0            M  (  1 )  =  0 

.c 

X: 

SHEAR 

NCMENT 

SLOPE 

DEFLECTION 

EXACT    DEF. 

0.00 

-0.124999 

C.COOOOO 

0. CI  49 72 

C. 000000 

0.000000 

0.05 

-0.  124999 

-0.0Co250 

0.014316 

C. 000746 

0.000746 

0.10 

-0.124999 

-0.01250  0 

0.C14347 

C. 00147b 

0.001477 

0.15 

-0.  12  499  9 

-C. 016750 

0.C13566 

C.002175 

0.002176 

0.20 

-0.124999 

-0.025000 

0 .012472 

C. 002 32 7 

0.002*28 

0.25 

-J. 124999 

-0.021250 

0.C11066 

C. 003417 

0.003418 

0.30 

-0. 119999 

-0.037415 

0  .C093-+3 

C. 003929 

0.0  039  20 

0.35 

-0.134999 

-Q.04:>030 

0.007333 

C. 004347 

0.00<t322 

0.40 

-0.07S999 

-C. 047745 

0.C05O57 

C. 004653 

0. J04609 

0.45 

-0.0*4999 

-0.050910 

0.002534 

C. 004349 

0  .004766 

0.50 

O.GJCOOO 

-0.052075 

0.000000 

C. 004914 

0. 00<t785 

0.55 

0.C45CJ0 

-C.0509i.Q 

-0.C02534 

C. 004-349 

0.0047o6 

O.oC 

0. C3CC00 

-0.047745 

-0  .005057 

C.G04o53 

0.0046C9 

O.o5 

0.105000 

-0.043030 

-O.C07333 

C. 004347 

0.004  322 

0.7C 

0.12CC00 

-C. 037415 

-0.009343 

C. 003929 

0.003920 

0.75 

0.125000 

-0.031250 

-0.C11066 

C. 003417 

C. 003418 

0.3C 

0. 125000 

-0.025300 

-0.C12472 

C. 00^327 

0.002328 

0.85 

0. 125000 

-0.013750 

-0  .013565 

C.  002173 

0.002176 

0.9G 

0.125000 

-0.012500 

-0.C14347 

C. 001476 

0.001477 

0.95 

0.125COO 

-C.CG6250 

-0.C14316 

C. 0007^6 

0.000  74o 

1.00 

0.125000 

C.GQOOOO 

-0  .014972 

-C. 000000 

0.000000 
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fairly  well  even  in  this  case.   Nevertheless,  the  error  should 
be  expected  at  the  "free"  (unknown)  B.C.'s  since  the  other 
extremes  are  fixed.   Here  we  make  use  of  the  fact  that  the 
loading  is  symmetrical  and  the  results  should  be  symmetrical 
too.   This  is  actually  the  case  in  the  computer  printout,  so 
we  can  conclude  that  the  results  are  indeed  accurate. 

7.  Example  4 :   Clamped-Clamped  Beam 

Refer  to  Fig.  5. Id.   This  is  a  statically  indeterminate 
beam  clamped  at  both  ends  and  loaded  with  a  totally  arbitrary 
discontinuous  load  known  only  as  a  table  of  values  at  discrete 
points  as  shown  on  the  next  page.   We  shall  determine  the  shear, 
moment,  slope  and  deflection  of  the  beam. 

8 .  Solution 

Before  we  proceed  with  this  example,  let  us  remind 
ourselves  that  we  are  not  restricted  to  handle  only  exact- 
integrable  forcing  functions.   We  can  deal  with  arbitrary 
functions  as  well,  and  this  example  intends  to  illustrate  the 
procedure.   In  this  problem,  the  computer  program  reads  in 
data  from  a  table  of  values  instead  of  obtaining  them  by 
evaluating  a  given  function.   With  this  data  the  program  per- 
forms trapezoidal  integration  in  the  usual  way.   The  forcing 
function  is  defined  by  three  other  functions  as: 


f  (x)   =0  0  <_   x  <  0.3 


f2(x)   =   (from  data  deck)    0.3  <  x  <  0.7 


f  (x)   =   0  0.7  <  x  <  1.0 
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OAT  5    CSCK     FCF 
01  STANCE 
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.  JC 

b. 

4J0GG0 
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200000 

7. 

o  JoGoO 

6. 

9J0000 

6  . 

300000 

4. 

1000  00 

4. 

000000 

4. 

600000 

■a 
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4. 

QOGOOG 

3. 

200000 

3. 

^00000 

3. 

200000 

3. 

40 00 00 

3. 

OOOOGO 

C. 

J  JOG  00 

G, 

000000 

C. 
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c« 

GOOjGO 

0, 
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G. 

GOOOJO 
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C. 
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Unfortunately,  no  exact  solution  has  been  obtained  for 
this  case  and  we  are  not  able  to  compare  the  computer  results 
with  the  exact.   Nevertheless  some  conclusions  can  be  drawn. 
The  B.C.'s  for  a  clamped-clamped  beam  are  as  follows: 


Shear: 

Y"» 
o 

— 

Unknown 

Y  "  • 
n 

— 

Unknown 

Moment : 

Y" 
o 

= 

Unknown 

Y" 
n 

= 

Unknown 

Slope: 

Y' 

o 

= 

0 

Y 
n 

= 

0 

Deflection 


=   0 


n 


=   0 


To  determine  the  initial  conditions,  we  use  the  known 
relationships  and  substituting  in  the  known  values,  we  obtain 


0   =   h2S,  +  N,h2Y,M  +  nhY" 
1     1    o        o 


0   = 


h3S   +  N  h3Y,M  +  N,h2Y' 


o   o 


'1   o 


This  is  a  pair  of  simultaneous  algebraic  equations 
with  two  unknowns.   Solving  for  the  initial  conditions  we 
get: 


Y" 
o 


n(N  S0  -  nS  i 
o  2 o_ 

nN,  -  N 
1     o 


(5.8) 


Y  •" 

O 


S   -  N,S0 
o 1  2 

nN,  -  N 
1     o 


(5.9) 
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For  this  example  the  values  supplied  by  the  computer 
are : 


Shear:  Y'"  =   -0.948411 

o 


Moment:         Y"   =    0.230507 

o 


The  computer  printout  is  shown  on  the  next  page.   Here 
we  have  another  column  to  show  the  arbitrary  loading.   Note 
that  we  are  using  100  intervals,  but  actually  every  fifth  is 
printed. 

Perhaps  the  only  way  to  analyze  these  results  would  be 
to  check  if  the  B.C.'s  have  been  met  or  not.   The  results  show 
that  they  have.   Another  clue  to  assure  some  accuracy  would  be 
to  look  at  the  maximum  deflection  and  the  minimum  slope.   The 
loading  was  intentionally  concentrated  on  the  middle  of  the 
beam  so  we  can  expect  the  highest  deflection  and  minimum  slope 
in  this  neighborhood  since  we  have  similar  B.C.'s  at  both 
ends.   We  can  say  that  this  too  checks.   Consequently  we  see 
that  the  results  are  likely  to  be  trusted. 
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EXAMPLE    <t 
CLANFEC-CLAMPfC    3 -"AN 


B  CMC  A  FY 

CGMCI7  ir,\s  : 

Y (C)=0.0 
Y  '  ( 0 )  =  0  . 

Y(1)=0.0 

J         V  (  1)=0.0 

X: 

l:^c 

it-:.',R 

PQri?.  IT 

3LDPI 

DEFLECTIQ! 

0.  JC 

o. ccccoo 

-C.94c411 

0.230507 

C. JGOOuO 

O.OOuOOO 

0.  05 

0.0  JCCOO 

-C.  948411 

0  .133036 

C.  J10340 

0.  J0J263 

0.10 

0. CCOC JO 

-C.S^d^lx 

0  .  135o66 

C.C13309 

J. JO 0  9 94 

0.15 

0. CCCCOO 

-C.9434J.1 

0.C38245 

C. 023  906 

0.J02059 

0.20 

o. j:co30 

-L.94c41 1 

0  .040825 

C.0^71j3 

0.003344 

0.25 

c. ^ccooo 

-C. 9484 11 

-O.C0fc596 

C. 027939 

0.00^732 

U.3C 

1.  L  '«.  u  G  o  J 

-c 943411 

-0.C53991 

C.02O474 

0.0061C3 

0.^5 

2.3ococo 

-v-  .  C^ri^l  1 

-O.C990o7 

C.  02262  3 

0.  J07339 

0.4C 

o. w^COjJ 

-0.6 7 jti i 

-0.1374a7 

C . Jl 6630 

0.J0833G 

0.45 

6.4:ccoo 

-0.35  ?4ll 

-0  .163 t28 

C. J09094 

0 . 306  979 

0.  50 

o  .  i  J  0  0  J  0 

~u»uJu91 i 

-0  .173116 

C.u006i4 

U. 00^224 

0.  2b 

6.7:ccoo 

0.279086 

-0  .16  7032 

-C.JO  79 52 

0.009  0  39 

U.oC 

6  •  J  vj  u  Ju  0 

C. 609085 

-0.  14  4  9  53 

-C. 01531  i 

0.  J08*40 

0.  6ia 

A-  .  3  c  ccoo 

o  .  c:255  J  j 

-0  .  1033  4 1 

-C. 0^2^07 

0.007^33 

0.  70 

3. CJCCoJ 

L  .99o^30 

-J  .C63133 

-C. J26542 

0.036255 

0.  75 

j • u jCLJU 

1.C1 153J 

-0  .012634 

-C.02344J 

0.004370 

0.3C 

J.  CC  CC JO 

x.Qli.530 

0  .  C  3  7  3  ^  > 

-C. J27310 

0.U03454 

0.d5 

J . C  J  cooo 

L. 011580 

0  .  C  8  8  4  72 

-C. 024651 

0.00  2  132 

G  .  9  C 

0 • 0  j  OOOO 

i. 01 1530 

0.139050 

-C.01 39o3 

0.001032 

J. 95 

0. ccccoo 

1  .01 i530 

0.139t>29 

-C. 010  746 

0.000279 

1.  00 

0.0JCC00 

1.011580 

0.240207 

Q.JOOOOO 

-C. 000000 
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THE    COMPUTER    PROGRAM 


file:    slope  »atfiv        ai 

tJCE  XREF 

C*** ********************************************** ********** 
C***************  THE  COMPUTER  PROGRAM  ******************** 
C* ************************************************ ********** 

c 

C  THIS  PROGRAM  CALCULATES  THE  SHEADING  FCRCE.  6ENCING  MOMENT.  SLOPE  AND 

C  DEFLECTION  OF  EEAMS  BY  PERFCRMNG  FCUR  SUCCESIVE  INTEGRATIONS  ON 

C  THE  LCAOING  FUNCTION.  THE  FCLUOWING  CONFIGURATIONS  ARE  CONSIDERED: 

C  A. -CANTILEVER  ( CLAMPEC-FRE E )  BEAM 

C  2.-CLAMPED-ROLLER  BEAM 

C  3. -SIMPLE  SUPCRTED  (FIN-ROLLER)  EEAH 

C  4.-CLAMPED-CLAMPED  BEAK 

C  THE  LOADING  CAN  3E  ANY  NUMBER  OF  0ISTRU6UTED  ANC/CR  CCNCENTKATEO 

C  FORCES.  CONCENTRATED  MOMENTS  CAN  BE  OECCMPCSED  INTO  A  CCUPLE  AT  THE 

C  MINIMUN  POSSIELE  DISTANCE  ACTING  AT  THE  PCINT  CF  APPLICATION  CF  THE 

C  MOMENT. 

C 

C  THE  PROGRAM  CONSISTS  OF  THREE  ROUTINES: 

C  1.  SUBROUTINE  SLOPE 

C  2.  SUBROUTINE  INCCM 

C  3.  SUBROUTINE  LCAC 

C  WHICH  ARE  DEFINED   IN  THE  FOLLOWING  PARAGRAPHS 


C 

C********  MAIN  PROGRAM  *************************** 

C 

C  READS  IN  NUMBER  CF  STEPS  DEFINED  BY  THE  INTEGER  VARIALBLE  N 

C  CALLS  SUBROUTINE  INCCM  FOR  INITIAL  CCNDITICNS 

C  CALLS  SUBROUTINE  SLCPE  FOR  FINAL  RESULTS 

C 

C********  SUBROUTINE  INCCM  *********************** 

C 

C  READS     IN     THE      INTEGER     VARIABLE    K     WHICH     SPECIFIES     THE     KIND     OF    PROBLEM 

C  WITH     THE     FOLLOWING     CODE: 

C  K=I     CANTILEVER     fclEAN     (INITIAL     CONDITIONS     MUST     BE     SPECIFIEC) 

C  K  =  2    CLAMPED     ROLLER     BEAM 

C  K=3    CLAMPED-CLAMPED     BEAM 

C  K=4    SIMPLE     SUPPCRTEC     BEAM 

C 

C     THE  PARAMENTERS  ARE  DEFINEO  AS  FOLLOWS: 

C  CALL  INC3N(N.  YC  .Y 1C .Y20 . Y3C J 

C     N  =  NUMBER  CF  STEPS 

C  YO    =  INITAL  CEFLECTICN 

C  YIO  =  INITIAL  SLOFE 

C  Y«C  =  INITIAL  BENDING  MOMENT 

C  Y30  =  INITIAL  SHERING  FCRCE 

C 

C  THIS  SUBROUTINE  SUPPLIES  THE  INITIAL  CONDITICNS  BY  CALLING  SLCPE. 

C  FCR  A  CANTILEVER  EEAM  THIS  STATEMENT  IS  NCT  EXECUTED 

C 

C*********  SUBROUTINE  SLOFE  ********************** 

C  PERFORM  TWO  TASKS:   1.  CALCULATES  AND  SUPPLY    INCCM  WITH  THE  NECESARY 

C  PARAMETERS   TC  DETERMINE  INITIAL  CCNCITICNS 

C  2.  EXECUTES  THE  SOLUTION  OF  THE  PROBLEM 

C 

C  THIS  SUdROUTINE   IS  A  THE  COOING  OF   THE  RELATIONSHIPS  GIVEN  BY 

C  EOLATIONS  3.13  TO  3.32.  THE  SUMMATION  TERMS  ARE  CONTAINED  IN 

C  THE  COMMON  STATEMENT 

C 

C  THE  CALLING  PARAMETERS  ARE  DEFINEO  AS  FULLOWS: 

C  CALL  SLCPE ( N. YO  .Y 1C ,Y23 . Y3C. J  J 

C     N  =  SAME  AS  IN  INCCM 

C  YC   =  " 

C  VIC  =  " 

C  Y2C  =  " 

C  V3C  =  " 

C     J  a  TAKE  TWO  VALUES:  J=0:   TC  FIND  INITIAL  CCNDITICNS  ONLY 

C  J=i:  TC  EXECUTE  THE  SOLUTION 

C 

C********  SUBROUTINE  LOAD  ************************ 

C 

C  SUPPLIES  SLOPE  WITH  THE  NECESARY   INFORMATION  GBTAINEO  FRCN  THE 

C  LOADING  FUNCTICN.   IT  CAN  EE  AN  ACTUAL  FUNCTION  CR  A  DATA  CECK. 
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file:  slope     »«tfiv    ai 

c  must  de  written  specifically  for  each  n£w  prc8le*.  when  dealing 
c  with  concentpatco  loads.  the  variable  rim  in  routine  slope  is 
c  is  set  equal  tc  the  concentrated  lcao  directly.  at  the  point  cf 

C  ACTICN. 

C  THE  CALLING  PARAMETERS  ARE  CEFINED  AS  FCLLCWS 

C  CALL  LOAO!e.H2) 

C  U-    CURRENT     VALLE     OF     THE     VARIABLE     X     1LENGHT)     GIVEN     eY     RCUTINE     SLOPE 

C  H2=    CONTAINS     THE     VALUE     CF     THE    LOADING     FUNCTION     AT     X=8 

C 

C*** *************************************** ******* 

C  MAIN     PROGRAM 

C***** ************************************ ******** 

REAL     YO. Y 1C. Y2C. Y30 

INTEGER     J.N 

COMMON    SO.S  t  .  S2.  SNO.SM  .RN.H 

READ!  5. I  0  )N 
10     FORMAT!  I  J  ) 

H=l ./FLJAT(N-I } 

PN=FLJAT( N-l ) 

CALL      INCON(N. YO. YlO. Y20 . Y30 ) 

J=0 

CALL     SLOPE  (N.  YO,  YlO.  Y  £0  ,  Y20.J  ) 

STOP 

END 
C 
c*** ******************************* ************** 

SUBROUTINE     LGA0(e.H2) 

REAL  a.ri2 
F (X  )  =  .  ..  . 
H2=F( 6) 

RETURN 
END 
C 

SUBROUTINE     INCCN(N.Y0.Y1C.Y20.Y3C) 

REAL     Y0. Y 1C. Y20.  Y30 
INTEGER     J.K.N 

COMMON     SO  .S  I  .  S2.  SNO.SM  .RN.H 
REAO!5.10)K,YC.Y10.Y2C.Y3C 
10     FORMAT)  I  1  ,4F  10. A  ) 

IF( (K.EQ. Cl.Ofi. IK.GT.4I)      GO     TC     60 
IF(  K.EO.  I  )     GC     TC     SO 
J  =  l 

CALL     SLGPElN. YC. YlO. Y20. Y30. J ) 
IF( K-3)20 .30.40 
C*****  *****************************  *************** 
C  CLAMPEO-RCLLER 

C***  *****************************  ***************** 
20    WRITE16.2E) 

25     FORMAT!  •  1  •. 'EXAMPLE     2  •  / / I  X .  'C LA MPEO- ROLL ER     EEA M' // 1  X .• BOUNDARY     CON 
*DITIONS:  V(0)=0.0  Y!l )=0.0'//27X .' Y      (0)=0.0  Mil)=0.0'//) 

Y2O=H*(S2*SNO-RN*S0 )/(RN*SNl-SNO) 
Y30-         (S0-SN1*S2) /( RN*SN l-SNO ) 
RETURN 

C  CLAMPEO-CLAMPEO 

30     WRI TE(6.35 ) 

35     FORMAT!  •  1  •.  'EXAMPLE     4  •  / / 1  X .  'C L A MPED-C L A MPED     BE  AM • // I  X .  • BOUND Afi Y     CO 
*NDITIONS:  Y(0)  =  0.0  Y(  1) =0.0* //25X.  'Y     (0)=0.0  Y     il)=0.0'//) 

Y20=H*( S0*SNl-Sl*SN0)/(fiN*SNO-SNl*SNl ) 

YJO=        ( S1*SN1-RN*S0 )/ (RN*SNO-SNl*SNl ) 

RETURN 

C  ^     PIN-ROLLEfi^ 

40     WRI TE16.4E ) 

45     FORMAT!  •  1  •,  'EXAMPLE     3 ' / / 1  X .  '°  I N-POLLEfi     B E A M • //  I X .  'BOUNDARY     CONOITI 
*ONS:  Y!0)=0.0  Y! 1 )=0.0 '//27x. 'M( 0)=0 .0  Mil)=0.0'//) 

Y10=(  (H/RN )**2)*( 3N0*S2-RN*S0  ) 
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file:    slope  »atfiv        ai 

Y30  =    -S^/fiN 
RETURN 

C  CL4MPED-FPEE 

5C     iHITE(6,5S)»0,yiC.V2a.r2C 

55    FORMAT< •  1  •  .  'EXAMPLE     1  ■ // I  X . "C LAMPEO-FREE     BEAK*         //  I  X . • BCUNO AR Y     CON 

COITIONS:  Y<  0) =•  .F8  .4.4X.  • Y     (0 )  =  • ,F8.4 .4X.//27X,      •*( 0 )  =  • .F 8 .4 . 

*4X. •V|0|=«.F8.4//| 
GO     TO     90 

60     *RITF(6.6EJ 

65     F09MAT|iX,'N     FRQK      1     TC     4     CNLY • / ) 
90     RETURN 
ENO 

c 

SUBROUTINE     SLCPE(N.Y0.Y10.Y20.Y30.J) 
C*** **************************************** ****** 

REAL     SU00.SU0I.SU02.SU03.SU11 . SU 12 . SU 1 3 . SU22 . SU23 .SU33 . RI N T . 
*RINTl.RINT2.RINT3.Hl.h2.RI.a.SNMl.SNM2 

INTEGER     I.J.K.L.N.M.NS1K .NSOL .NSIL.IM1 

COMMON     SO  .SI  .52.  SNO.SM  .RN.H 

SU0O=0. 

SUO 1=0. 

SU0  2  =  0. 

SU03=0. 

SU1 1 =0. 

SU1 2=0. 

SUl 3=0. 

SU22=0. 

SU23=0. 

SU33=0. 

R  IMT  =  0. 

RINT1 =0. 

RINT2=0. 

RIMT3=0. 

NS1 K=0 

NS0L=0 

NS1L=0 

H1=0.0 

M2=0.0 

IF(  J.GT.O  )     GO     TO     100 

WRl TE(6.2C) 
100     00     400     1=1. N 

IM1 =1-1 

RI=FLOAT(   IM1  ) 

B=RI*H 

CALL     UQA018.H2) 
IF»   I  .EG.  1  )GC     TO     200 

RINT  =  0.5*I-*<H1+H2J 
200     RINT4=RINT3 

RINT3=R[NT2 

RINT2  =  RI NT  1 

RINT1=RI NT 
C 

SU33=SU33+RINT1 

S3        =SU33 
C 

SU23=SU23*RINT2 

SU22=SU22*SU23 

S2=.5*SU33+SU22 
C 

<  =  I-2 

IF(  I.LE. 2  )     K  =  0 

SUl 3=SUl 3+RINT3 
SUl 2=SU12+SL 13 

SUl 1=SU1 I *SU12 

NS1 K=NS1 K *K 

SNM l=FLQ AT(NS1K ) 

SN1      =0.5*R I ♦SMMl 

S1  =  .25*SU3  3*SL22«-SU1  1 
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file:    slope  *atfiv         ai 


L  =  I-3 

IF(  I.LE.  3 )     L  =  0 
SU03=SUQ3*R INT4 
SU0  2=SUC2*SU0  3 
SU0  1=SUO 1 +  SU02 
SUOO=SUGO  ^SCIO  1 
NSIL=NSIL+L 
NSOL=NSOL*NSlL 
SNM AFLOAT  (NSOL  ) 
SNO     =0.25SR H-SNM I +SNM2 

SO=. 1 2S*SU33*.7S*SU22+1 ,5<SU1 H-SUOO 
IF( J.GT.O )     GO     TO     300 
C 

Y3=S3>V3Q 

Y2=H* ( S2  +  R I*Y30) »V2C 

Yl  =  l(SH-SM?Y3C)*HtfiI*Y2C  )*W-Y10 
Y  =  <   (  (S0*SN0*Yja|*H  +  SNl*Y2C)'H  +  BI*YI0)*H+YC 
*RITE(6,  tC»b.V3.Y2.tt  .  * 
300     MI=M2 
400     CONTINUE 

IF(J.GT.O)     GO     TO     999 
WRI  TE(6. 99  ) 
10     FORMAT!  1 X  ,F4  .2 .2X.4F 1 2.6/  ) 

2  0    FORMATI 2X .«X: ».     gx.'SHEAB'.bX. 'MCMENT '.6X. • SLOPE* .5X. *OEF LECTION* ) 
99    FORMAT)  •  I  •  ) 
999     RETURN 
END 
C 
SENTRY 
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VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

Based  upon  the  research  carried  out  in  this  thesis  and 
the  results  obtained,  the  following  conclusions  can  be  drawn: 

1)  First  of  all,  the  method,  in  the  way  it  has  been 
developed,  shows  the  fact  that  the  solutions  are  totally 
independent  of  each  other.   Any  value  of  the  unknown  can  be 
found  without  generating  all  previous  ones.   A  closer  look  at 
equations  (3.13)  to  (3.16)  reveals  that  we  can  apply  any  of 
these  relationships  at  any  point  i  (0  <  i  <  n)  directly,  pro- 
vided we  know  the  initial  conditions  and  the  summation  terms 
which  can  be  generated  in  advance.   In  all  of  these  equations 
there  is  only  one  summation  term  which  depends  on  the  given 
function  f(x),  the  other  summations  are  series  of  integers 
totally  independent  of  the  given  problem.   This  fact  represents 
a  good  saving  in  terms  of  computer  time  if  it  is  conveniently 
exploited . 

2)  As  it  has  been  shown  in  the  examples,  the  power  of 
the  method  is  perhaps  its  ability  to  deal  with  arbitrary  func- 
tions and  this  is  important  since  many  engineering  problems 
lead  to  these  kinds  of  functions. 

3)  So  far,  there  has  not  been  any  error  correction  in  a 
strict  sense.   The  only  corrective  measure  has  been  step  size 
reduction.   However,  equation  (4.9)  shows  that  correction  can 
be  performed  provided  we  know  how  the  local  error  behaves. 
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When  integrating  constant  and/or  linear  functions,  (4.9) 
applies,  but  for  other  cases  it  does  not.   In  any  case,  equation 
(4.3)  indicates  that  the  local  error  depends  only  on  the  second 
integral  which  in  turn  is  one  degree  higher  than  the  given 
function  itself.   This  suggests  the  idea  that  if  we  know  how 
the  original  function  behaves,  linearly,  quadratically ,  etc., 
we  could,  to  some  extent,  predict  the  error  behavior  and  carry 
out  a  correction.   Further  research  is  recommended  in  this 
particular  case. 

Nevertheless,  as  shown  in  the  examples,  some  accuracy  has 
been  achieved  even  working  in  single  precision.   For  more 
complicated  problems  involving  complex  functions  and  requir- 
ing very  accurate  solutions,  double  precision  is  still  a  good 
possibility. 
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